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Quantum electrodynamics (QED), the gauged quantum field theory of the electromag- 
netic interaction, has celebrated ground-breaking successes in the twentieth century. In 
fact, its multi-loop predictions for the anomalous magnetic moments of the electron and 
the muon were found to agree with highest-precision measurements within a few parts of 
1CT 12 and 10 -10 , respectively. 

Another ultrapure laboratory for high-precision tests of QED is provided by positron- 
ium (Ps) , the lightest known atom, being the electromagnetic bound state of the electron 
e~ and the positron e + , which was discovered in the year 1951 pQ. In fact, thanks to the 
smallness of the electron mass m relative to typical hadronic mass scales, its theoretical 
description is not plagued by strong-interaction uncertainties and its properties, such as 
decay widths and energy levels, can be calculated perturbatively in non-relativistic QED 
(NRQED) [2], as expansions in Sommerfeld's fine-structure constant a, with very high 
precision. 

Ps comes in two ground states, 1 Sq parapositronium (p-Ps) and 3 Si orthopositronium 
(o-Ps), which decay to two and three photons, respectively. In this Letter, we are con- 
cerned with the lifetime of o-Ps, which has been the subject of a vast number of theoretical 
and experimental investigations. Its first measurement [3] was performed later in the year 
1951 and agreed well with its lowest-order (LO) prediction of 1949 [I]. Its first precision 
measurement [5], of 1968, had to wait nine years to be compared with the first complete 
one-loop calculation [6] , which came two decades after the analogous calculation for p-Ps 
[7j being considerably simpler owing to the two-body final state. In the year 1987, the 
Ann Arbor group [8] published a measurement that exceeded the best theoretical predic- 
tion available then by more than ten experimental standard deviations. This so-called 
o-Ps lifetime puzzle triggered an avalanche of both experimental and theoretical activities, 
which eventually resulted in what now appears to be the resolution of this puzzle. In fact, 
the 2003 measurements at Ann Arbor [Hj and Tokyo [10] , 



r(Ann Arbor) = 7.0404(10 stat.)(8 syst.) /is" 
r(Tokyo) = 7.0396(12 stat.)(ll syst.) fis 
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agree mutually and with the present theoretical prediction, 



r(theory) = 7.039979(11) /is 
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The latter is evaluated from 




(3) 



where [4] 
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is the LO result. The leading logarithmically enhanced 0(a 2 \na) and 0(a 3 In 2 a) terms 
were found in Refs. [TT][T2] and Ref. [13], respectively. The coefficients A = —10.286606(10) 
[6lfTT1fT4lfT5lfT6j . B = 45.06(26) [15], and C = -5.51702455(23) pi] are only available in 
numerical form so far. Comprehensive reviews of the present experimental and theoretical 
status of o-Ps may be found in Ref. [T8] . 

Given the fundamental importance of Ps for atomic and particle physics, it is desir- 
able to complete our knowledge of the QED prediction in Eq. ([3]). Since the theoretical 
uncertainty is presently dominated by the errors in the numerical evaluations of the co- 
efficients A, B, and C, it is an urgent task to find them in analytical form, in terms of 
transcendental numbers, which can be evaluated with arbitrary precision. In this Letter, 
this is achieved for A and C. The case of B is beyond the scope of presently available 
technology, since it involves two-loop five-point functions to be integrated over a three- 
body phase space. The quest for an analytic expression for A is a topic of old vintage: 
about 25 years ago, some of the simpler contributions to A, due to self-energy and outer 
and inner vertex corrections, were obtained analytically [19], but further progress then 
soon came to a grinding halt. The sustained endeavor of the community to improve the 
numerical accuracy of A [6|ll|14|15p6] is now finally brought to a termination. 




Figure 1: Feynman diagrams contributing to the total decay width of o-Ps at 0(a). Self- 
energy diagrams are not shown. Dashed and solid lines represent photons and electrons, 
respectively. 

The 0(a) contribution in Eq. (JBJ), Ti = ToAa/ir, is due to the Feynman diagrams 
where a virtual photon is attached in all possible ways to the tree-level diagrams, with 
three real photons linked to an open electron line, and the electron box diagrams with an 
e + e~ annihilation vertex connected to one of the photons being virtual (see Fig. [I]). Taking 
the interference with the tree-level diagrams, imposing e + e~ threshold kinematics, and 
performing the loop and angular integrations, one obtains the two-dimensional integral 
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representation [T6_ 



ma 
36^2 



dxi dx 2 dx 3 



X\ x 2 x 3 



6(2 — x 1 — x 2 — x 3 ) 



x [F(xi, x 3 ) + perm. 



(5) 



where Xi, with < Xi < 1, is the energy of photon i in the o-Ps rest frame normalized by 
its maximum value, the delta function ensures energy conservation, and perm, stands for 
the other five permutations of Xi,x 2 ,x 3 . The function F(xi,x 3 ) is given by 



F(x 1 ,x 3 ) = g (xi,x 3 ) + 9i(xi, x 3 )hi(x x , x 3 ), 



(6) 



i=l 



where gi are ratios of polynomials, which are listed in Eqs. (A5a)-(A5h) of Ref . |I6] , and 
/ii(xi) = ln(2xi), h 2 {x 1 ) = JzrOi, 



h 3 (x 1 ) 
hi{xi) 
K(x x ,x 3 ) 
h 7 (x h x 3 ) 



-L[C 2 -Li 2 (l-2x 1 )] ) 

^[3C 2 -2^], hix^^-el, 

, = _ [2 Li 2 (r+, 00 - 2 Li 2 (r B , ^ 

Z^f X\X\X 3 X 3 1 

- Li 2 (r+,0) + Li 2 (r^,0) 



with Xi = 1 — x,- and 



#i = arctan(yxi/xi), #i = arctan(yxi/xi), 

r\ = yJW[ 1 ± 



IXiX 3 



X\X 3 



B 



i- 



Here, (2 = vr 2 /6 and 



(7) 



(8) 



Li 2 (v 



~J — ln(l - 2rt cos 6» + r 2 t 2 ^ 



(9) 



is the real part of the dilogarithm [see line below Eq. (|20|) ] of complex argument z = re l9 
[2D] . Since we are dealing here with a single-scale problem, Eq. ([5]) yields just one number. 

Although Bose symmetry is manifest in Eq. ([5]), its evaluation is complicated by 
the fact that, for a given order of integration, individual permutations yield divergent 
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integrals, which have to cancel in their combination. In order to avoid such a proliferation 
of terms, we introduce a regularization parameter, 5, in such a way that the symmetry 
unter Xi Xj for any pair i ^ j is retained. In this way, Eq. (jSJ) collapses to 

7 1-5 1-5 

„ ma f , /" dx 2 s ,„„ s 

OTT^ J J £1X9X3 

25 l-xj+5 

where X3 = 2 — X! — x 2 - Note that we may now exploit the freedom to choose any pair of 
variables Xi and Xj (i 7^ j) as the arguments of F and as the integration variables. 

The analytical integration of Eq. (fTUj) is rather tedious and requires a number of tricks 
to be conceived of. For lack of space, we can only outline here a few examples. Specifi- 
cally, we consider the last two functions of Eq. ([7]), which are most complicated. Using 
Eq. and after some manipulations, we obtain the following integral representation for 
h 7 (xi,x 3 ): 



1 



4 J y/t(x!X 3 - XiX 3 t) 

+ 2 ln(x3 + x 3 t) — lilt 



In — — 

X1X3 

(ID 



Exploiting the X\ «-> X3 symmetry of the coefficient g 7 (xi,X3) multiplying hr(xi, x 3 ), this 
can be simplified as 

-In*]. (12) 

At this point, it is useful to change the order of integrations. Observing that the log- 
arithmic terms in Eq. (1121) are Xi independent, we first integrate over X\ (for a similar 
approach, see Ref. [2T]). In order to avoid the appearance of complicated functions in the 
intermediate results, the integration over t in Eq. ffl2l is performed last. 
Analogously, h 6 (xi,x 3 ) can be rewritten as 

*.(*„*,) = -£ /-^ -[In^-lnz, 

2 J y/t{XiX 3 - X\X 3 t) 

+ hi(x 3 + x 3 t)}, (13) 

in which the part proportional to lnxi and the complementary part are first integrated 
over x 3 and x±, respectively. The t integration is again performed last. 

Let us now consider a typical integral that arises upon the first integration: 

1 1 

I = Jj J ^ ln[l - 4t(l - t)(l - a;)] ln(l - x). (14) 
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Direct integration over t or x would lead to rather complicated functions in the remaining 
variable. Instead, we Taylor expand the first logarithm using ln(l — x) = — Y^=i x n /n to 
obtain 



00 4 n r dt r dr 

/ = - E - / y^ 1 - *r J ^(1 - *) n hi - *). 

n—1 n n 



(15) 



Now the two integrals are separated and can be solved in terms of Euler's Gamma function, 
T(x) = f™dte-H x -\ Using 



dx 



x 



(l-x) n ln(l -x) = -ip'(n + 1), 



(16) 



where if)(x) = dlnr(x)/dx is the digamma function, we finally have 

I 



00 4" T 2 (n) 



iy 2n r(2n) 



(17) 



Another class of typical integrals yields sums involving digamma functions of half- 
integer arguments, e.g. 







ln[l + 4t(l -t)(l- x)] ln(l - x) 



x-2 



-A) n T 2 (n) 



E 



'0' 



, f n + 2 



1/ 



, fn+ 1 



/ and J belong to the class of so-called inverse central binomial sums [221123] , and methods 
for their summation are elaborated in Ref. [23j. With their help, / and J can be expressed 
in terms of known irrational constants, as 



-4C 2 / 2 2 -|-8Li 4 
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-c 



4, 



3 2 /| 2 2 11 3 47 4 

J — ~2^ 2 + "4 ~ ^C 2 ^r + + y^^r + 288^ 

+ 4Z 2 Z r Li 2 (r) + ^ 2 Li 2 (r) - 6l 2 Li 3 (-r) 
o 

4 (\ 
- 2/ r Li 3 (-r) + 5Z 2 Li 3 (r) + -/ r Li 3 (r) + 6 Li 4 - 



+ 4Li 4 (-r) - 5Li 4 (r) 



y^Si, 2 (r) + ^Si, 2 (r 2 
19 



- 4 S 2 , 2 (-r) + 5 S 2 , 2 (r) + (3/2 + — 

o 

where r = (y/2 - l)/(v / 2 + 1), k = lnz, 



Sn,p(^') 



(n — 1)! p! 7o t 



(19) 



(20) 



is the generalized poly-logarithm, \A n {x) = S n -i,i(x) is the poly-logarithm of order n, and 
Cn — C{ n ) — Li ra (l), with ((x) being Riemann's zeta function [201125] . 

Unfortunately, not all integrals can be computed so straightforwardly. In more com- 
plicated cases, the integrations are not separated after expansion into infinite series. We 
then rely on the PSLQ algorithm [24] . which allows one to reconstruct the representation 
of a numerical result known to very high precision in terms of a linear combination of a 
set of irrational constants with rational coefficients, if that set is known beforehand. The 
experience gained with the explicit solution of the simpler integrals helps us to exhaust 
the relevant sets. In order for PSLQ to work in our applications, the numerical values of 
the integrals must be known up to typically 150 decimal figures. 

After a laborious calculation, we obtain 

2, o rtN , 56 901 , 11303 , 19, 2701 , , 

-(tt 2 -9)A= C2 C4 + —k (2^2 

9 V ; 27 216 192 s 6 108 

253 251 913 83 21 

+ "24 C2 * 2 + 144 2 + IP 2 * 3 + 256 3 " T Wr 

21 7 63 

- y^2 Li 3 (-r) - -Z r Li 3 (-r) + —l 2 U 3 (r) 

63, T . . . 249 T . / 1\ 249 T . fl 
+ ^ Ll3(r) -^2- Ll4 (-3) + T6- L H3 

+ ^ Li 4 Q + 7Li 4 (-r) - 7S 2 , 2 (-r) 

63 t . / \ 63 . , 11449^ 91 , 

- -j L W) + — S 2 ,2(r) + -^~Cs - yCsfe 



35^ , 1 



yC2^-^/ 3 -y/rLi 2 (r) 



35 t . / \ 175 . . 14 , 2 , 119 

+ — Li 3 r " Sl > 2 r + ~T Si, 2 r 2 + — ( 
6 3 3 3 

The constant C in Eq. ([3]) is related to A through [17] 



(21) 



A 229 

C = - - — + 8l 2 . (22) 
3 30 V ; 

From Eqs. (!2~T]) and (J22l) . A and C can be numerically evaluated with arbitrary preci- 
sion, 

A = -10.28661 48086 28262 24015 01692 10991 ... , 

C = -5.51702 74917 29858 27137 88660 98665 .... (23) 

These numbers agree with the best existing numerical evaluations [T6lfl5] within the 
quoted errors. 

In conclusion, we obtained the 0(a) and 0(a 3 In a) corrections to the total decay 
width of o-Ps, i.e. the coefficients A and C in Eq. (jBJ), respectively, in closed analytic 
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form. Another important result is the appearance of new irrational constants in Eq. (I2"T1) . 
These constants enlarge the class of the known constants in single-scale problems. The 
constant B in Eq. ([3]) still remains analytically unknown. 

We are grateful to G.S. Adkins for providing us with the computer code employed for 
the numerical analysis in Ref. [16] . This work was supported in part by BMBF Grant 
No. 05 HT6GUA, DFG Grant No. SFB 676, and HGF Grant No. NG-VH-008. 
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